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ABSTRACT 

Using the full, three-dimensional potential of galaxy cluster halos (drawn 
from an N-body simulation of the current, most favored cosmology), the dis- 
tribution of the X-ray emitting gas is found by assuming a polytropic equation 



of state and hydrostatic equilibrium, with constraints from conservation of en- 
ergy and pressure balance at the cluster boundary. The resulting properties of 
the gas for these simulated redshift zero clusters (the temperature distribution, 
mass-temperature and luminosity-temperature relations, and the gas fraction) 
are compared with observations in the X-ray of nearby clusters. The observed 
properties are reproduced only under the assumption that substantial energy in- 
jection from non-gravitational sources has occurred. Our model does not specify 
the source, but star formation and AGN may be capable of providing this energy, 
which amounts to 3-5 xl0~ 5 of the rest mass in stars (assuming ten percent of 
the gas initially in the cluster forms stars). With the method described here it is 
possible to generate realistic X-ray and Sunyaev-Zel'dovich cluster maps and cat- 
alogs from N-body simulations, with the distributions of internal halo properties 
(and their trends with mass, location, and time) taken into account. 



Subject headings: cosmology:theory — galaxies:clusters:general — intergalactic 
medium — X-rays:galaxies:clusters 
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Introduction 



Current and upcoming surveys in a variety of wavelength bands will increase the num- 
ber of well-observed clusters of galaxies by at least an order of magnitude, while probing 
to much higher redshifts than before. Understanding the physical state of the intra-cluster 
medium (ICM) will be essential to exploiting this new data. In particular, it is necessary to 
develop methods of accurately modeling the thermal state of the gas in clusters before one 
can extract cosmological information from large surveys, which measure quantities arising 
from that state. For cluster-sized halos in a cosmological setting, the theoretical final distri- 
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However, the hot intr- 
acluster gas in these systems does not parallel the DM in either density or temperature 
distribution. 

Much progress has been made in understanding the expected ICM distri bution inside a 



standard DM halo (with the density profile showing a pow er law cusp as in iNavarro et al. 



19971 ; iMoore et al.lll999l . or similar). iMakino et al.l (119981 ) gave an analytic expression for 
the density of isothermal gas in hydrostatic equilibrium with a NFW p otential; this was 



1998 



soon extended to non-isothermal gas with a polytropic equation of state (ISuto et al 
Wu et al.ll2~0~00l; iLoewensteinl |2000| ; lAscasibar et al.l 120031 ). and to triaxial halos (ILee &; Sutol 
20031 : IWang fc FanlboOfj ). The resulting gas profiles possess a finite density core, and not a 
cusp as seen in the DM. 

The gross energetics of the gas do not parallel that of the DM either. Assuming that 
the gas energy comes solely from gravitational collapse gives the self-s imilar scalings between 



Eke et al. 



mass M, luminosity L, and temperature T of M oc T 3//2 and L cx T 2 ( IKaiserlll986 
19981) . However, these scalings do not agree with the observed relations, leading [Kaiser 



19911 ) to propose that non-gravitational energy injection is important. This idea has gained 
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ICM in DM potentials 
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gas with s hort cooling times, low entropy, and low total energy, leaving behind higher entropy 
mate rial (jVoit &: Bryan! l200ll ; iTozzi &: Norman! l200ll ; IVoit et al.l 120021 ; IScannapieco fc Oh 
2004 ). Detailed computer simulations including star for mation and feedback are confirm- 
ing the importance of non-gravitational processes (e.g. iBorgani et al.l 120051: lEttori et al.l 
20061 : ISiiacki & SpringeJbooafBorgani et al.lbood : iRomeo et alihood : iMuanwong et alJliooS 



Nagail 120061 . and references therein). It appears from the cited papers that both processes 
are at work: low entropy gas is incorporated into stars, and energy and metals are added to 
the remaining gas via feedback processes. 

Based on these advances in our understanding of clusters, one can usefully combine 
prescriptions for gas physics with analytic modeling of DM profiles, but this method has 
limitations. Since halos are formed b y stochastic merging of s ubunits, there is true scatter in 
halo concentration and inner slope (Avila- Reese et al. 19991; Jin^2000 ; Subramanian et al. 



2000; Bullock et al. 2001 
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2006 
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2006|, and references therein) and substructure (e.g. lGill et al.ll2004l ; lReed et al 
200d . and references therein). The clustering of halos h as been shown 



to depend upon age or concentration, as well as mass (jBerlind et ajj|2006j , and references 
therein). Thus mass, environment, and formation history all play significant roles in deter- 
mining halo properties. 



In prior work, lOstriker et al.l (120051 ) improved on the use of analytic DM potentials by 
instead using the full three-dimensional potential of halos drawn from N-body simulations, 
combining these detailed 3-D models with current modeling of the gas physics. This pro- 
cedure has the advantage of including the full distribution of halo concentration, as well as 
halo triaxiality and substructure. By drawing on a large N-body simulation volume (compu- 
tationally much less costly than a full hydrodynamic run), trends of internal halo properties 
with mass, location, and time are all included, along with halo-halo correlations. This pro- 
cedure inevitably requires the utilization of certain dimensionless parameters which — since 
they derive from feedback processes — are difficult, at present, to determine from ab initio 
computations. The purpose of the present paper is to tie down th ese parameters using ob- 
servations of X-ray clusters. In particular, we apply the method of lOstriker et al.l (120051 ) to 
a set of halo s drawn from a larg e DM simulation of a cosmology in accord with the WMAP 
3-year data (jSpergel et al.ll2006l ). The resulting catalog is compared to X-ray observations 
of nearby clusters. The amount of energy input from non-gravitational sources can signifi- 
cantly affect gas properties, but with the proper amount — consistent with AGN activity — 
this method can reproduce the properties of the local cluster population. 
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The procedures used to create the simulated cluster sample are described in $2], these 
clusters are compared to X-ray observations in §3], and we discuss the implications in §|4j One 
extensi on of this wo r k is t o use the gas properties to compute the Sunyaev-Zel'dovich (SZ) 
effect; ISehgal et al.l (120061 ) do this to make available simulated large- area, sub-arcminute 
resolution microwave sky maps. 



2. The Simulated Cluster Catalog 
2.1. Dark Matter Halos 



To produce a population of DM halos, we chose cosmological parameters to match 
the results from the 3-year WMAP data combined with large-scale structure observations 
( jSpergel et al.l 120061 ) . The spatially flat LCDM model was used, with total matter density 
Q m = 0.26, baryon density Qb = 0.044 (so the cosmic mean baryon fraction f c m 0.169). 
and cosmological constant Qa = 0.74. Also, the Hubble constant Ho = 72kms _1 Mpc~ 1 
(i.e. h = 0.72 = H /100kms _1 Mpc _1 ), the primordial scalar spectral index n s = 0.95, and 
the linear matter power spectrum amplitude erg = 0.77. Concerning the N-body simulation 
parameters, the number of particles N = 1024 3 and the box size L = 1000/i~ 1 Mpc, making 
the particle mass m p = 6.72 x lO lo /i- 1 M ; the cubic spline fsee iHernquist fc Katzl \198% 
softening length e = 16 .3/t~ 1 kpc. The init ial conditions for the iV-body run were created with 
the GRAFIC1 code^ ( Bertschinger 2001 ). with a few modifications. Because the spherical 
Hanning filter employed in this code to isotropize small-scale structure also significantly 
suppresses power on small scales, it w as not used. The linear DM transfer function at z = 
was calculated with the CMBFAST JZaldarriaga fc Seliakl boOdh codfl The DM growth 
factor was used to scale the resulting power spectrum to the initial simulation redshift, 
chosen to be when the density fluctuation amplitude on the grid scale is 10%, z = 35.3. 
With these parameters, a cluster with mass ~ 7 x 10 14 /i _1 M Q would contain 10 4 particles; 
a typical core radius for such a cluster would be ~ 250/i _1 kpc, or 15 times the particle 
softening length. 



The simulation was run with the TPM cod^l flBode et al.ll2000l ; Bode fc Ostrikerl 120031 ). 
with a couple of improvements over the publicly released version. Most variables in the 
improved code are double precision (including particle positions and velocities), with the 



1 This code is available at http://web.mit.edu/edbert/ 
2 Available at http : //www . cmbf ast . org/ 



3 Available at 



http : //www . astro .princeton . edu/$\sim$bode/TPM/ 
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main exception of accelerations and potentials, which are still single precision. Also, no 
lower limit was set to the p arameter B used in the TPM domain decomposition (see Eqn. 



5 of IBode fe Ostrikerl 120031 ). This means that at late times there will be more particles 
followed at full force resolution, leading to improved simulation of the lowest mass objects; 
at z = all cells with eight or more particles were followed with trees. The initial domain 
decomposition parameters were A = 1.9 and B = 9.2, the PM mesh contained 2048 3 cells, 
and the maximum sub-box size was 256 cells. At the end of the run, 54% of the particles 
contained in 2% of the total volume were being followed with 5 x 10 6 trees. 

The standard friends-of-friends (FOF) halo finder with linking length b = 0.2 was run 
on the simulation volume at z — 0, identifying almost 2 x 10 6 halos with 30 or more par- 



ticles. The resulting mass function agrees well with the fitting formula of IWarren et al. 



( 120061 ); the difference in the cumulative mass function is less than 5% over the range 
2 x 10 12 < M/ o/ /(/i -1 M ) < 10 14 and less than 15% above this. In what follows, only 
halos containing gas with temperatures above 1.5keV will be included in our discussion. We 
find only halos with M500 > 1O 13 /i _1 M will ever meet this limit (though typically it takes 
over twice this mass). 1O 13 /i _1 M corresponds to 150 particles, so the question arises as 
to whether or not there is sufficient numerical resolution for such halos. To test this, we 
compared the halos described in this paper with a sample taken from a higher resolution 
N-body simulation. This high resolution run has a box size L = 320/i~ 1 Mpc and a spline 
softening length e = 3.2/i _1 kpc, so the mass resolution is thus improved by a factor of 30 
and the spatial resolution by a factor of 5; otherwise it was generated and evolved in exactly 
the same manner as the simulation described above. When applying the gas prescription of 
the next section, we used a mesh cell size of I = 9.60/i -1 kpc, a factor of 3.4 smaller than the 
standard case (it was difficult to go to any finer mesh because of the memory requirements of 
the resulting large computational arrays). The maximum feedback model (described below) 
was assumed. Tabled] compares the mean and standard deviation for several observables pre- 
dicted in the standard and high resolution runs, using two mass bins corresponding roughly 
to temperatures of 1.7 and 2 keV; little difference is evident. A smaller cell size leads to 
a smaller volume around each halo (for the same number of cells); this accounts for most 
of the differences seen in the two samples — the mass resolution is relatively unimportant. 
Increasing the cell size to / = 12.80/i _1 kpc in the high resolution run results in a distribution 
even closer to the standard run. 
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2.2. The Gas Prescription 



The gas distribution in each halo is calculated according to the prescription of lOstriker et al 



( 120051 ). Gas is placed in hydrostatic equilibrium with the DM gravitational potential of the 
halo using a polytropic equation of state. Pressure balance with infalling gas near the virial 
radius and energy conservation determine the two constants required for the polytropic fit. 
Two important processes alter the gas energy. Star formation removes low entropy gas; we 
fix the conversion of gas into stars at 10%. This leaves the most important free parameter, 
which is the energy input into the cluster gas via feedback processes. We will show that with 
a reasonable amount of feedback it is possible to match X-ray observations of hot cluster 
gas. 

In detail, a cubic mesh enclosing the particles is placed around the halo, with the cell 
size twice the particle spline softening length, or I = 32.55/i -1 kpc. The mass m p for each 
particle is placed on the mesh using the cloud-in-cell method, yielding the DM density, 
Pd/cj for each cell k. The gravitational potential on the mesh, is computed from the 
density as in a standard Particle-Mesh code, but with a nonperiodic FFT. The center of the 
cluster is defined as the cell with the lowest potential, 0o = MIN(<fik). The radii enclosing 
various overdensities are calculated at this point, the outermost being the virial radius, r vir , 
enclosing the overdensity expected from spherical tophat collapse, or 97 times the critical 



density at z = for the cosmology used here (this is a change from lOstriker et al.l 12005 
where overdensity 200 was used). The velocity of the cluster as a whole is taken to be the 
mean velocity of the 125 particles closest to the cluster center (or, for halos with fewer than 
250 particles, the innermost half). Particle velocities are moved to the rest frame of the 
cluster, and then the kinetic energy (KE) of each particle is placed on the grid in the same 
manner as the mass, yielding the KE per unit volume |iDfc- 

It is assumed that gas originally had the same distribution as the DM, with density 
fcPvk and KE f c \tr>k (fc = Qb/^m)- A certain amount of the gas mass, M* (described 
below), will have turned into stars; this is presumably the most bound material, so cells 
are ranked by binding energy cj) k + |tofc, and then cells are checked off until the sum of the 
masses f c pr>kl 3 equals M*. The initial mass M g and energy E g of the remaining gas are thus: 

M g = J2fcPDkl 3 , (1) 

k 

Eg = J2fA ( t ) kPDk + lt Dk }l 3 , (2) 

k 

where the sum is over all cells inside r V i r except those marked off for star formation. Also, 
the gas surface pressure P s on the cluster exerted by surrounding material is estimated from 
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the kinetic energy in a buffer region nine cells (293/i 1 kpc) thick outside of r % 



iv r 1 Yl ^ fi>L 



k=l 



fc. 

3 



(3) 



where the sum is over the Nb cells in the buffer region r v i r < < r v i r + 91. 



Now suppose the gas is allowed to rearrange itself within the DM potential such that 
it is in hydrostatic equilibrium and has a polytropic equation of state with index T = 1.2. 
There is much support for such a model — see the discussion and Fig. 1 in lOstriker et al. 
(120051 . comparing the p olytropic model with a full, high resolution cosmological simulation 
by G. Bryan), and also lAscasibar et al.l (120061 ) . We will treat the gas as a tracer such that 
the potential, set by the DM, does not change. Defining 



0* 



1 



r 



Pa 



[l + 5 rel )TP 



(4) 



as m 



Ostriker et al.l (120051 ). the resulting gas pressure P and density p are given by 



P{fk) 



Poor 

- Poof 1 



(5a) 
(5b) 



Here 5 re i is a nonthermal component of pressure, assumed to be proportional to thermal 
pressure such that the total P to t = (1 + 5 re i)P. To specify the final gas distribution given 
these assumptions, two quantities still need to be determined, namely the pressure P and 
density po at the potential minimum. This can be done with two equations of constraint, 
derived by requiring conservation of energy and by matching the external surface pressure, 
as follows. For a given choice of Po an d Po, the final radius rj of the gas initially inside 
be found by summing outwards from the cluster center until the initial mass M g is 
enclosed: 

£ PoO^l 3 = M g . (6) 

r k <rf 

This implies that gas may expand or contract, changing the gas fraction inside r vir . Assuming 
the external surface pressure changes little with radius, there will be mechanical work done, 
causing a change in energy proportional to the change in volume, AE P = (47r/3)(r^ r — rj)P s . 
The equation for conservation of energy is thus 



rk<rf 



i 

r-i 
k ' 



+ ^(i + 25 rel )p ef' 



l 3 = E g + AE P + e f M*c 2 



(7) 
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The term e/M*c 2 is feedback inferred from supernovae and AGN, discussed in more detail 
below. Matching the final surface pressure to the external pressure yields the other equation 
of constraint: 



k=l 



again summing in a buffer region Tf < r < r^ + 9/. (Note the S re i terms were omitted in 
Ostriker et al.ll2005l ). With Eqns. (GO) and (jHJ) it is possible to iterate (e.g. with Newton- 
Raphson) to a solution for the final gas density and pressure (or temperature) profile. 

For this paper we will assume that at z = the initial (that is, inside r vir prior to any 
rearrangement) star to gas ratio is 10%, in other words /* = M*/M g = 0.10, which implies 
M+ = ( f r M vi r ) f*/(\ + /*)• This ratio agrees well with the value in n earby clusters measured 



bv iLin et al.l (120031 ). and is slightly lower than that measured by IVoevodkin fc Vikhlinin 
( 120041 ) . Star and black hole formation will return energy to the remaining gas via super- 
novae and AGN activity — writing this energy as tfM^c 2 of course assumes this energy is 
proportional to stellar mass. This seems plausible: the number of supernovae is expected 
to be proportional to the mass in stars, and the mass of cent ral black holes in spheri- 



ods is roughly proportional t o the stellar mass in these systems ( iMerritt &: Ferraresd 12001 
Kormendy &: Gebhardtll200ll ). Im proved observat i onal constraints may alter this assump- 



tion. The rough estimate given in lOstriker et al.l (120051 ) is that ej « 3 x 10~ 6 for SN and 
€f fa 4 x 10~ 5 for AGN, so we will take ej = 5 x 10~ 5 as the maximum case. This is roughly 3 
keV per particle for the gas inside the virial radius, which is at the high end of the plausible 
range. 



2.3. Cluster Temperature 



To characterize the temperature of the gas we will use T ew , the X-ray emission- weighted 
T inside a projec ted radius of Rmn - The X-ray luminosity is calculated using the cooling 
function A(T) of iMaller fc Bullock! J2004h for T < 10 8 A", and assuming A(T) cx T°- 5 fo r 



T > l O 8 -?^; the metallicity is set to one third solar (IBaumgartner et al.ll2005l ). iMazzotta et al. 



(120041 ) showed that the projected spectroscopic temperature of a thermally complex cluster 
will in fact be lower than the emission-weighed temperature. However, this difference will be 
more pronounced when computing a single-temperature fit to a full hydrodynamic simulation 
(containing shocks, cold fronts, and other short-lived structures) than it is in our simple 
equilibrium models (which lack such local inhomogeneities) because such features, which 



increase the thermal comple xity of the gas, contribute to this systematic bias (IMazzotta et al. 



20041 : iKawahara et al.ll2007l ). In order to quantify this effect, we compared T ew with T sp , the 
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spectroscopic temperature mea sured in the range 0.15i?soo < R < -Rsoo- We used the code^ 
developed by IVikhlininl (120061 ) to compute T sp , using the Chandra response function and 
Galactic absorption Njj = 2 x 10 20 cm -2 . Note that simply excluding the center reduces the 
measured T sp relative to T ew , independent of spectral effects. For the maximum feedback 
model, we find kT sp = 0.94kT ew — 0.06keV, with little scatter; in other words, the two agree 
within 10%. Given this small difference, we will use the conceptually simpler T ew unless 
stated otherwise. 



3. Effects of Feedback 



3.1. Gas Temperature 



X-ray surveys provide valuable information on the luminosity, temperature, and mass 
of clusters. In this section we explore these properties as derived using the method of the 
previous section. The relationship between mass and temperature is shown in Fig.[TJ The 
temperature is the emission- weighted T ew inside a projected radius -R50O) an d the mass is 
that contained in a spherical overdensity of 500 times critical, r 50Q . Lines show the median 
mass at a given temperature, and the shaded regions enclose 68% of the clusters; the cases 
of no feedback (e/ = 0) and maximum feedback (e/ = 5 x 10~ 5 ) are shown. Feedback has 
little effect for the hottest, most massive clusters: the feedback energy is small compared 
to the binding energy of these clusters, and thus of little importance to the dynamical state 
of the gas. Feedback has a greater effect in less massive clusters, making the gas somewhat 
hotter. Still, at M 500 w 5 x lO^h^M®, feedback increases T by only ~ 33%. One effect 
not apparent from the figure is that for masses below M500 ~3x 10 13 /i~ 1 M Q the maximum 
feedback can be enough to make the total gas energy positive, unbinding the gas from the 
halo. Thus our method produces no halos with temperatures below abou t 1.5 keV in the 



maxim um feedback case. The data points shown in this figure are from McCarthy et al. 



(200 4j), who combined ASCA observations with the extended ROSAT HIFLUGCS sample of 
Reiprich &: Bohringerl (120021 ); cluster cores were not excluded when T was determined. Only 
those clusters closer than z < 0.06 are shown. Below 4 keV, the feedback model provides 
a superior fit; both models are in agreement with the data above this. However, there is 
significantly more scatter in the observed M — T relation than is produced in our model. 
Cooling (which we neglect beyond that involved in star fo rmation) will increase the scatter 
in this relation ([McCarthy et al.ll2004l ; lO'Hara et al.ll2006l ). The existence of young systems 
which are out of dynamical equilibrium can also broaden the observed M—T relation, but the 



4 Available at 



http : //hea- www. harvard. edu/$\sim$alexey/mixT 
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effects may not be very pronounced (jO'Hara et al.ll2006l ); we are to some extent accounting 
for this, because merging halos will have greater kinetic energy per particle and thus a higher 
temperature than a relaxed halo of the same mass. We have not modeled observational error, 
which will also of course add to any intrinsic scatter. 

Our predicted M — T relations change little if we instead use the spectroscopic tem- 
perature, as show n in Fig. El Also sh own are ten nearby, relaxed clusters observed with 
XMM-Newton by lArnaud et al.l (120051 ) and ten relaxed clusters observed with Chandra by 
Vikhlinin et al.1 (120061 ). Both used spectroscopic temperatures, and excluded the cores (al- 
though the radial range used to determine T is slightly different). These observations exhibit 
considerably smaller scatter, lending credence to the idea that differing dynamical states and 
cooling in cores will increase the scatter in temperature at a given mass. As was the case 
before, both models agree reasonably well with the observed M — T relation above 4 keV, 
but the high feedback case is a better fit in the 2-3 keV range. The M — T relation can be 
well fit by the power law 



E(z) 



Mr 



500 



A 



© 



kr y 

5keVy 



(9) 



where E(z) = H(z ) /Ha. lArnaud et all (120051) found A = 2.69 ±0.10 and a = 1.71 ±0.09 for 
their sample, while I Vikhlinin et al.l J2006h obtained A = 2.89±0.15, a = 1.58±0.11. We first 
attempted to fit this relation to the simulated halos with ordinary least-squares regression 
in the log-log plane, but this was dominated by the more numerous low mass halos and 
unduly affected by outliers. Thus we instead adopted the following procedure: we divided 
the x-axis into 20 logarithmically spaced bins, calculated the median in each bin, and then 
found the best fit to these points. This in effect gives more massive clusters a higher weight 
in the fitting; the resulting fits follow closely the median lines in the figures. As listed in 
Table [2], fitting to all halos with kT sp > 3keV in the zero feedback model gives A = 3.12 
and a = 1.49 ± 0.02. This follows the self-similar slope of 1.5, and gives cooler clusters at 
a fixed mass than is observed. In the maximum feedback case temperatures shift to higher 
values, and the slope becomes steeper: A = 2.56 and a = 1.62 ± 0.03. This slope agrees well 
with the observations, although the normalization yields slightly hotter clusters at a given 
mass. The formal error we obtain for A is small (near one percent) so in Table [2] we give the 
rms fractional difference of the halos from the best-fit relation; this scatter reflects well the 
width of the shaded regions in the figures. If the amount of feedback €f varied from cluster 
to cluster, then the scatter seen would be larger. 

Another method of characterizing nearby clusters is the temperature function, which 
does not require a mass determination. Th e distribution of cluste r temperatures is sensitively 
dependent on the cosmological model; in lOstriker et al.l (120051 ). which used the first-year 
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WMAP power spectrum amplitude, the fit to observations was inadequate. Shown in Fig. [3] 
are two measurements of the cumul ative temperature function, with different methods of 
determining the cluster t emperature: llkebe et al.l (120021 ) excluded cluster cores when fitting 
for T, while iHenryl (120041 ) did not. For purposes of comparison, we took from the simulation 
a "light cone" covering one octant of the sky out to z = 0.2. This covers the redshift range 
used in the observations, although they are not volume limited. During the simulation, the 
matter distribution in a series of thin shells was saved; the radius of each shell corresponds 
to the light travel time from a z = observer sitting at the origin of the box, and its width 
corresponds to the time interval between shells. Thus a volume-limited mass distribution, 
including time evolution, is obtained. Locating halos and adding gas was done in the same 
manner as before. To compute the star/gas ratio at z > 0, the star format i on ra te was 
assumed to follow a delayed exponential model (Eqn. 1 of iNagamine et al.l 120061 ). with 
decay time r = lGyr. Both T ew and T sp are shown for the ej = 5 x 10~ 5 feedback model 
in Fig. [31 Because our cosmological parameters were chosen in part to match large-scale 
structure observations, and the simulated M500 — T relation matches that of nearby clusters, 
it is not surprising that our simulated temperature function is a reasonable fit to the observed 
one . Our XL,, wh ich includes the core, gives a higher abundance in the 3 — 6 keV range than 
the iHenrvl (12004 ) data , while our T sp , excluding the core, instead gives a lower abundance 
than llkebe et al.l (120021 ) . The zero feedback model appears to give temperatures that are too 
low; thus based simply on T, it appears that some non-gravitational energy input is required 
to explain the properties of existing clusters. The model with feedback and WMAP 3-year 
cosmological parameters now provides a good fit to the observed temperature function. 



3.2. Gas Density 

Other cluster observables are more dependent on the gas density, most notably X-ray 
luminosity. The top panel of Fig. H] shows the L x —T relation of our simulated catalog for three 
values of feedback, again with medians shown as lines and shaded regions enclosing 68% of 
the clusters; here L x is the bolometric luminosity inside a projected radius of -Rsoo- Unlike the 
M — T relation, feedback produces a significant change in L x at a given temperature because 
the Bremsstrahlung emission is proportional to the square of the gas density, but otherwise 
there are similar trends seen. Again the effect of feedback is less important in the most 
massive clusters, where gravitational binding energy dominates. Also, the scatter for a given 
amount of feed back is much smaller than that observed in nearby clusters (the data points 



are again from iMcCarthy et al.l 12004] . with z < 0.06). The zero feedback and e/ = 5 x 10~ 5 
models bracket the range of luminosities seen in nearby clusters. An intermediate model 
with €f = 3 x 10~ 5 is also shown in Fig.Hl it appears this is an insufficient amount to explain 
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the lowest luminosity clusters. Fitting a power law L x /(10 h ergs 



A(kT sp /5keV) c 



to all our halos with kT sp > 3keV yields a = 1.66 ± 0.03 in the zero feedback case, not as 
steep as the self-similar expectations of of L x oc T 2 . Including feedback ste epens this relation 



considerably to a = 2.51 ±0.08, more in line with th e observed value (e.g. lArnaud fc Evrard 



19991 ; iReiprich fc Bohringerl 120021 ; llkebe et al.ll2002l ). The exact value of the slope we find 
depends on the lower temperature limit used. 

A more direct probe of gas density is the gas fraction within a given radius. It appears 
that the gas fraction increases with increasing radius, and higher temper a ture clusters have 



a higher gas fraction as well (IDavid et al.l Il995l ; lArnaud &: Evrardl Il999l ; iMohr et al.l 11999 
Vikhlinin et al.l 120061 ) . The gas fraction from our model is shown in the bottom panel of 



Fig.lU with f g defined as the fraction of the total mass inside a spherical radius r 500 enclosing 
an overdensity 500 times critical; for the purpose of computing the total mass, we assumed 
that stellar mass followed the same radial profile as the dark matter. The model curves 
display the type of behavior one might expect based on the top panel: models with feedback 
show significantly lower gas fractions (i.e. lower densities at a given T), with the effect being 
most pronounced for the lowest mass clusters. Recent observations of t he gas fraction by 



Vikhlinin et al.l (120061 ) are shown in the figure, as well as two clusters from lGastaldello et al. 



wh ere for these po i nts th e temperature was derived from M500 using the M — T 
relation of IVikhlinin et al.l (120061 . note they use a spectroscopic temperature, and exclude 
the central regio n when find i ng T) . Also shown as a dashed line is the best fit f g — T 
relation found by lMohr et al.l (119991 ) using 45 ROSAT clusters and temperatures taken from 
the literature. Here again, some feedback is required to bring the models in agreement with 
observed gas fractions, but again the spread in any given model is too small to fit all observed 
values of f g . 

Also shown as a dotted line in the figure is the mean gas fraction from our cosmological 
model, after turning enough gas into stars to make the global star/gas ratio ten percent, or 
fg = fc'(l — f*/ (! + /*))• Without any feedback, the baryon fraction inside r 50 o will reflect 
the cosmic mean value, but energy input drives this fraction lower, particularly for smaller 
clusters. We find the gas fraction increases with radius, so at overdensities higher than 500 
this discrepancy will be even greater. This raises the question of how far out one must go 
before the cluster contains a fair sample of the cosmic mass budget. Fig. shows the median 
gas fraction inside the virial radius as a function of temperature (still using T ew inside -R500); 
the virial overdensity is 97 times critical for our chosen cosmology. Even at this radius, it 
is only for the most massive clusters (kT ew > 6keV) that the baryon fraction reaches the 
cosmic mean; in smaller clusters feedback causes the gas to expand, reducing the gas fraction 
by many tens of percent. Attempts to determine the cosmic baryon density from clusters 
will need to take this effect into account. 
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Fig.[5]also demonstrates how the median gas fraction changes in the maximum feedback 
case if we also add a relativistic component with 5 re i = 0.20 (there is little change in the 
spread around the median). With this component, a lower temperature is required to achieve 
pressure balance at a given density. This model behaves like the no feedback case at higher 
masses, and like an intermediate feedback case at lower masses. This behavior also holds 
for all the other relationships (M 50 o — T ew ,L x — T ew , etc.) explored in this paper; thus it 
seems a relativistic component will have little effect on thermal cluster observables. While 
our implement ation is quite simplified , similar results were found using full hydrodynamic 



simulations by iPfrommer et al.l (120061 ): they found that including the effects of cosmic rays 
caused only small changes in the gas fraction and integrated SZ signal (they found a larger 
change in L x , but this was related to cooling cores, which we do not implement). 

To summarize this section, we find that a WM AP 3-year cosmological model coupled 
with a feedback parameter of e/ = 3 — 5 x 10~ 5 (which corresponds to an input energy of 
roughly 2-3 keV per baryon) provides a good fit to the extant X-ray observations of hot gas 
in clusters. 



4. Discussion 

In this paper we have presented a method for determining the gas distribution inside a 
fully three-dimensional potential; this method assumes hydrostatic equilibrium and a poly- 
tropic equation of state, and also that the original gas energy per unit mass equals that 
of the DM. We then applied this method to a z = catalog of DM cluster halos drawn 
from N-body simulations, and compared the resulting ICM distributions to observations of 
nearby clusters. The main result is that this simple gas prescription can reproduce many 
of the observed bulk properties of the ICM, including the temperature distribution and the 
relationships between temperature and mass, X-ray luminosity, or gas fraction. The main 
drawback is that in nearby X-ray clusters there is significantly more scatter seen in these 
relationships than is produced in our method. This could be due to a number of factors, 
including cooling cores, our assumption of hydrostatic equilibrium, and observational errors. 

The advantages of using this type of model for the ICM are clear. It is possible to 
simulate a large volume with a N-body code at much smaller computational cost than is 
required for a full hydrodynamical treatment. When adding gas with the method described 
here, the distributions of internal halo properties — concentration, triaxiality, substructure, 
etc. — are taken into account, as are their trends with mass, location, and time, plus any 
alignments and correlations between halos. However, it should be noted that this approach 
does not account for the dynamical effects of a baryonic component. Results from halo 
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formation simulations demonstrate that including a dissipational gas co mponent can alter 
the radial profile of the DM halo (IGnedin et al.l l2004j ; iLin et al.l 120061 ) and its ellipticity 
( IKazantzidis et al.l 120041 ). which would in turn alter the ICM distribution. 



A second result is that a significant amount of non-gravitational energy input is re- 



the same conclusion (e.g. 


Kaiser 


1991: Balosm et al. 


1999: i 


3uto et al. 19981: Wu et al. 


20001: 


Locwcnstcin 


2000: Tozzi & NormankoOl: Komatsu k Seliak 


2001 


:lBabul et al. 2002: Voit et al. 


20021: Dos Santos & Dore 


2002; 


Shimizu et al. 2004: 


Lapi et al. 


20051: lAfshordi et al. 


2005: 


Solanes et al 


20051). In this paper we do not trv to find the sources of this input. Instead. 



we simply ascertain what level of feedback would produce the observed relations among 
cluster observables. Processes involved in star formation provide some of this energy, but 
not enough. The most likely source for extra energy is AGN, which can conceivably deliver 
enough feedback to explain the temperature, X-ray luminosity, and gas fraction distribu- 
tions of local clusters. However, there is little margin for error: if we have overestimated the 
amount of star formation, the black hole to stellar mass ratio is less than 0.0013, and/or the 
amount of energy returned to the gas by black hole formation is less than 3% of the black 
hole rest mass, then the required energy of 2 to 3 kev per pa r ticle will not be produced. 
This conclusion differs somewhat from that of iMcCarthy et al.l (120061 ). who find that AGN 
heating is an implausible (though not impossible) explanation of cluster gas fractions; they 
calculate that in order to reduce the gas fraction within rsoo to the observed level it would 
take 10 keV per particle in the gas observed within r 50 o (roughly 0.12 M 5Q0 ). Rescaling 
this value to our normalization (energy per particle of the gas mass initially inside r vir , or 
f c M V i r ), makes this a required energy of roughly 4 keV per particle. The reason we require 
less energy is in part due to a different initial state: examining the lower panel of Fig.HJ by 
just accounting for star formation but not including any feedback, one can see our method 



l eads to a gas fraction inside r 50Q already lower than the cosmic mean, while McCarthy et al. 



(120061 ) started with a state in which the gas fraction equals the cosmic mean. (Although in 



any case they argue that efficiencies of AGN outbursts are ~ 10 3 , not ~ 1 x ). Running hy- 



drodynamic simulations without radiation or feedback, ICrain et al.l (120061 ) found the baryon 
fraction inside r^n of 90%, which agrees well with our method at higher masses. However, 
Crain et al.l (120061 ) also find this fraction still holds at lower masses, and that the baryon 
fraction is still 90% at the virial radius (again with no feedback), while we find a higher 
fraction in both of these instances. Note that if we had instead started with an initial state 
consistent with a 90% baryon fraction, then we would require a lower amount of feedback to 
reproduce observed gas fractions. 



It is useful in this context to compare with recent hydro dynamical simulations which 
include feedback. The left panel of Fig. [6] relates gas to total mass inside r 50 o; the solid line 
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shows the best- f it pow er law relation at z = found from adaptive refinement simulations by 
Kravtsov et al.l (120061 ). With no feedback we find a similar slope (slightly larger than unity), 
but for a given halo mass there is a higher gas fraction than in the simulations. Adding 
feedback reduces this discrepancy at higher ma sses, but at lower m asses our maximum 
feedback case instead has lower gas fractions. As iKravtsov et al.l (120061 ) did not attempt to 
include AGN feedback it is not surprising that an intermediate amount of feedback would 
provide the best match. We do not match the simulations whe n it co mes to the total mass- 



temperature relation, howe yer. For Eqn. ([H I), IKravtsov et al.l (120061 ) find A = 3.85 ± 0.19 
and a = 1.524 ± 0.07, while iKav elTall (bood i with an SPH simulation find A = 4.47 ± 0.19 
and a = 1.76 ± 0.07. We obtain a similar slope (a = 1.62), but more importantly we find a 
different normalization — that is, we find hotter temperatures at a given mass, as would be 
expected because we included a higher level of feedback, 

The near future will see a number of surveys that select clusters of galaxies via their 
Sunyaev-Zel'dovich (SZ) decrement, which is proportional to the gas pressure i n the cluster 
integ rated along the line of sight. Currently the Sunyaev -Zel'dovich Array (S ZA; ICarlstrom et al 
2000l ) and the Arc Minute Microkelvin Imager (AMI; iKneissl et al.l l200ll ) are equipped to 
perform such a search on tens of square degrees on the sky. However, the Atacama Cosmol- 
ogy Telescope (ACT), the South Pole Telescope (SPT), APEX, and ulti mately the Planck 



surveyor wil l scan thousands of square degrees on the sky in the radio (IRuhl et al.l 12004 



Fowlerl 120041 ; iGiisten et al.l 120061 : IClavel fc Tauberl 120051 ) . These surveys will detect thou- 
sands of clusters; for example, the SPT will scan 4,000 deg 2 on the sky and observe of the 
order of 6, 000 clusters of galaxies (this is for a flux sensitivity of about 1.5mJy at the 4-sigma 
detection threshold with a 1 arcmin beam operating at 150 GHz, and assumes as = 0.75). 
This translates to a limi ting mass of My m 10 14 ' 2 /i _1 M Q if one assumes the clusters are in 



hydrostatic equilibrium. ISehgal et al.l (120061 ) found a similar limit for a 90% complete cluster 
sample from ACT. 

The redshift distribution of clusters is very sensitive t o the amplitude and growth of lin- 



ear perturbations, and hence to cosmological parameters ( 


Holder et al. 


2001; 


Haiman et al. 


2001; 


Weller et al. 


2002; 


Battve & Weller 


2003; 


Maiumdar & Mohrll2004 


Youneer et al. 


2006) 



However, in order to exploit SZ cluster number counts one is required to understand the se- 
lection function of these surveys. This is most easily accomplished in terms of the flux 
decrem ent, which depends on the system temperature, exposure time, band-width, and effi- 
ciency (IBattye fc Weller! 120051 ) . In order to obtain cosmological constraints, it is necessary 
to convert the observables into a limiting mass of the survey. There are two approaches to 
obtain this mass limit. One is to start with the assumption that all clusters are spherical 
and in hydrostatic equilibriu m, and then include some nuisance parameters to allow for dev i- 
ations from this assumption (jVerde et al.ll2002l ; IBattye fc Wellerll2003l ; lYounger et al.ll2006l ). 
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Another approach is to use a very general parameterization of the mass-observabl e relation 



which in its most general form could easily introduce forty unknown parameters (IHul 12003 



Lima fc Hull2005l ). Currently there is little data to constrain the free parameters in either 
approach. In future o ne can exploit the SZ cluster observations themselves, to self-calibrate 
these free parameters (jHull2003l ; iLima fc Hull2005l ; iMajumdar fc Mohrll2004j ; iBattve fc Weller 
20031 ). However if one employs the most gene ral parameterization, little power is left in the 
surveys to constrain cosmological parameters (jHull2003l ; ILima fc Hull2005l ). Another possibil- 
ity would be to use compl ementary observations, such as weak lensing to cross-calibrate the 
mass - observable relation (IMajumdar fc Mohrll2004l ; iThe Dark Energy Survey Collaboration 
20051 ; ISealfon et al.l 120061 ) . A useful approach would be to have a physical parameterization 
of the mass-observable relation with some prior pro bability on the free p arameters and then 
self-calibrate the SZ surveys for these parameters ( Younger et al.l 120061 ) . However in order 
to obtain this prior knowledge we can not yet rely on observations because currently they 
are sparse, and in the near future observations will not resolve clusters because the beams 
of the instruments are typically larger than 1 arcmin. We hence require simulations to ex- 
plore the scatter in the mass-observable relation; in order to obtain realistic results, a large 
representative sample of galaxy clusters is required. 



Sehgal et al.l (120061 ) have already applied the method of this paper to a full light-cone 
N-body output (out to z = 3) in order to generate and make publicly available large-area, 
sub-arcminute resolution microwave sky maps. We intend to provide a detailed analysis 
of the mass-observable relation in a forthcoming paper, and will give here only a rough 
qualitative discussion. In particular, we can use our z = simulated catalog to explore 
how the amount of thermal energy in the gas will affect the SZ signal. One can express the 
strength of the integrated SZ flux as 



"SZ 



dA / pkTdl 



(10) 



where the integration limits are along the line of sight through the entire cluster, and over 
area out to projected radius R^oo- The right panel of Fig. [7] displays how the SZ signal varies 
with cluster mass in our model, for the zero and maximum feedback cases. At higher masses 
the relation has little dependence on feedback. Feedback reduces the SZ signal somewhat 
(as it results in gas being pushed out of the higher pressure cluster cores) and makes the 
relation steeper. Fitting to halos with M 500 > 10 14 /i _1 M Q , we find L$z oc M^q 2 with zero 
feedback. This is clo se to the self-similar slope of 5/3 predicted for spherical profiles (e.g. 
Reid fc Spergell 120061 ); apparently triaxiality and substructure have little effect on this re- 
lation. The slope steepens only sligh tly to 1.69 for ey = 5 x 10 -5 ; this is in reasonable 



agreemen t with the analyti c results oflReid fc Spergel ( 2006) and with hydrod ynamic sim 



ulations (IWhite et al.l 120021 ; Ida Silva et al. 



2004 



Motl et al. 



20051 : iNagail 120061 ). The exact 



-17- 



slope found for this relation depends on the lower mass limit used; if we included lower mass 
halos (or weighted the high mass halos less) the resulting slope would be steeper. The left 
panel of Fig. [7J shows how L S z varies with temperature. The relation is quite tight, and 
again steepens with increasing feedback; fitting to halos with kT sp > 3keV, L sz oc T s 2 p 43 with 



no feedback, the exponent increasing to 2. 77 for 



U- 



- sp 

5 x 10~ 5 . 



The zero feedback slope 
agrees well with the adiabatic simulation of iNagail (120061 ). but the feedback model is steeper, 
because we are putting in more energy. Increasing 6f from zero to e/ = 5 x 10~ 5 lowers the 
SZ signal at kT P „, = 5keV by 35%, which is similar to but slightly less than the effect seen 
by INagail (120061 ) between his adiabatic and star formation runs. A more detailed comparison 
is difficult because we are using the projected SZ signal; also there are differences in the 
cosmological parameters. 

This work makes it clear that allowance for feedback will be necessary if one is to 
utilize the upcoming SZ surveys for precision measurements of cosmological parameters. 
Fortunately, X-ray observations allow us to calibrate the feedback parameter; adequate fits 
to the data can be obtained if the ratio of energy input to stellar mass ise/ = 3 — 5x 10 -5 . 
Uncertainties in this parameter will propagate into uncertainties in the mass-flux decrement 
relation for SZ surveys. However, it can be seen in Fig. [7] that even the extreme case of 
reducing e/ to zero hardly changes this relation for clusters with masses above 2 x 10 14 h~ l M & . 
Moreover, it is not the scatter in the mass-observable relation which makes it difficult for 
future SZ surveys to constrain cosmolo gical parameters, but rather it is the uncertainty in 
the scatter which is the main problem (ILima fc Hull2005l ). We will use the methods of this 
paper to explore these issues in future work. 
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Table 1. Test of Resolution for low- mass Halos 



1 < M 500 < 4 4 < Af 500 < 6 

Standard High res. Standard High res. 



M500 


3.18 ±0.55 


3.29 ± 0.51 


4.85 ± 0.57 


4.81 ±0.56 


kTeu; 


1.68 ±0.19 


1.75 ± 0.21 


1.99 ± 0.24 


1.97 ±0.21 


L x 


1.14 ±0.76 


1.21 ± 1.02 


3.63 ± 1.99 


3.79 ±2.77 


Lsz 


3.31 ± 1.95 


3.16 ± 1.95 


7.13 ± 3.91 


6.24 ±2.94 



Note. — Mean and one standard deviation. Units of M500, kT t 
L x , L sz : 10 l3 h- 1 M e , kcV, 10 42 fe- 2 crgs- 1 , lO 12 /i- 1 keVM 



Table 2. Power Law Fit Parameters 









= 


e/ = 5x 10" 6 


Relation 


A 


a 


A 


a 




- M 500 


0.46 ± 9% 


0.67 ± 0.01 


0.55 ± 10% 


0.60 ± 0.01 


M g - 


M500 


1.44 ± 10% 


0.96 ± 0.01 


0.80 ± 12% 


1.17 ± 0.02 


L x — 


M500 


1.15 ± 26% 


1.13 ± 0.02 


0.31 ± 25% 


1.59 ± 0.03 


Lsz - 


- M 500 


0.43 ± 25% 


1.62 ± 0.02 


0.35 ± 31% 


1.69 ± 0.03 


M 500 




3.12 ± 14% 


1.49 ± 0.02 


2.56 ± 19% 


1.62 ± 0.03 


M g - 


" kT'sp 


4.37 ± 17% 


1.40 ± 0.03 


2.42 ± 24% 


1.84 ± 0.06 


L x — 




4.09 ± 29% 


1.66 ± 0.03 


1.36 ± 34% 


2.51 ± 0.08 


Lsz - 


— kTsp 


2.86 ± 15% 


2.43 ± 0.03 


1.85 ± 16% 


2.77 ± 0.04 



Note. — For each relation Y — X, the best fit parameters of the form 
Y/Y = A(X/X ) a arc given. For Y = (M 500 , kT sp , M g , L x , L sz ), Y = 
(W^h^MQ, 5keV, lO 1 ^- 1 ./^, 10 44 fe- 2 crg s" 1 , lO^h-^eV M Q ). When 
X = M500, the fit is to all halos with X > X = 1O 14 /i _1 M ; when X = kT sp , 
Xq = 5keV and the fit is to all halos with kT sp > 3keV. The scatter given for 
A is the rms fractional difference of the halos from the best-fit relation. 
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Fig. 1. — M — T ew relation for two values of the feedback parameter. T ew is the emission- 
weighted temperature from all material inside projected -Rsoo- Lines show the median 
va lue, and shaded region s enclose 68% of the clusters. Filled circles are data described 



in 



McCarthy et al.l (12004 ). using only z < 0.06 clusters. 
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Fig. 2. — M — T sp relation for two values of the feedback parameter. T sp is the spectroscopic 
temperature excluding the inner 0.15i?soo- Lines sh ow the median value, a nd shaded regions 
enclose 68% o f the c lusters. Filled circles are from IVikhlinin et al.l (120061 ) and crosses from 
Arnaud et aD J2005h . 
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Fig. 3. — X-ray temperature functi on. Cr o sses a re from llkebe et al.l (120021 ) (who excluded 
cluster cores) , and circles are from iHenryl (120001 ) (who kept them) . Lines are the volume 
limited z < 0.2 temperature function from simulated clusters: solid: emission-weighted 
T ew from all material inside projected -R500; dotted: spectroscopic T spec excluding the inner 
0.15i?5oo; dashed: emission-weighted T ew with no feedback. 
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Fig. 4. — Top panel: L x —T relation for three values of the feedback parameter. For the 
simulated z = clusters, the lines are the me dian and the shaded re gions enclose 68% of 
the clusters. Filled circles are data described in McCarthy et al.l (120041 ). using only z < 0.06 
clusters. Bottom pan el: gas fraction inside r^no- Points with error bars are data from 



Vikhlinin et al.l (120061 ) and iGastaldello et al.l (120061 ). and the dashed line is the best fit to 



45 ROSAT clusters by iMohr et al.l (119991 ) . The dotted line is the cosmic mean adjusted to 
make the global star/gas ratio 10%. 
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Fig. 5. — Gas fraction inside the virial radius r vir , normalized to the cosmic mean value 
(adjusted for star formation), as a function of temperature. Lines are the median and shaded 
regions enclose 68% of the clusters. The line without shading is the median for ej = 5 x 1CT 5 
and 5 re [ = 0.20. 
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Fig. 6. — Gas mass inside r 500 as a function of M 500 (left) and kTU ( right) . The straight 
lines are the full hydrodynamic simulation result from iKravtsov et al.l (120061 ). Lines are the 
median and shaded regions enclose 68% of the clusters. 
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Fig. 7. — Integrated SZ decrement, L sz = f dA f pkTdl, as a function of temperature (left) 
and mass (right), for two values of the feedback parameter. Lines show the median value, 
and shaded regions enclose 68% of the clusters. 



